Characterizing steady state and transient properties of reaction-diffusion systems 
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In the past the study of reaction-diffusion systems has greatly contributed to our understanding 
of the behavior of many-body systems far from equilibrium. In this paper we aim at characterizing 
the properties of diffusion limited reactions both in their steady states and out of stationarity. Many 
reaction-diffusion systems have the peculiarity that microscopic reversibility is broken such that their 
transient behavior can not be investigated through the study of most of the observables discussed 
in the literature. For this reason we analyze the transient properties of reaction-diffusion systems 
through a specific work observable that remains well defined even in the absence of microscopic 
reversibility and that obeys an exact detailed fluctuation relation in cases where detailed balance 
is fulfilled. We thereby drive the systems out of their nonequilibrium steady states through time- 
dependent reaction rates. Using a numerical exact method and computer simulations, we analyze 
fluctuation ratios of the probability distributions obtained during the forward and reversed processes. 
We show that the underlying microscopic dynamics gives rise to peculiarities in the configuration 
space trajectories, thereby yielding prominent features in the fluctuation ratios. 

PACS numbers: 05.40.-a,05.70.Ln,05.20.-y 



I. INTRODUCTION 

Understanding general properties of systems that are 
far from equilibrium remains one of the most challeng- 
ing problems in contemporary physics. In recent years 
some remarkable progress has been made in the study 
of fluctuations in nonequilibrium small systems, see, e.g. , 

^ajJi^MJS Um, tu m m m m, es m, 

YL&l , [l9l . [20 . |2U l22l l23l . |24| . This progress is mainly due 
to the formulation of various exact fluctuation and work 
theorems that provide very generic statements that hold 
true for large classes of systems. Some of the fluctuation 
theorems deal with the rate of entropy production, either 
in systems that are in a nonequilibrium steady state 0, Q 
or in systems that initially are in equilibrium before be- 
ing driven out of equilibrium by an external force 1, [1] . 
One also distinguishes between detailed 0, [ll], [14 1 and 
integral [J, [ljj [l5[ fluctuation theorems in cases where 
the system is driven from one steady state to another 
in finite time. In addition, work theorems [1, 9 re- 
late the free energy difference between two equilibrium 
states to the amount of work done during the switch- 
ing from one state to the other. Many of these theo- 
rems have been verified in various experimental settings 
0, [H, [H [H, [H HI HI Hi, thus illustrating their 
usefulness for characterizing nonequilibrium systems. 

In the past diffusion-limited reaction systems have 
been proven to be extremely useful in order to understand 
the generic behavior of many-body systems far from equi- 
librium. Especially, the study of these systems is at the 
origin of our under standingof the properties of nonequi- 
librium phase transitions [25L HH . Whereas phase tran- 
sitions in reaction-diffusion systems are by now very well 
characterized, this is quite different away from these spe- 
cial points. In the present work we discuss different ways 
of characterizing the steady state and transient proper- 
ties of generic reaction-diffusion systems, thus contribut- 



ing to a more comprehensive understanding of these im- 
portant systems. 

Recently discussed extensions of exact fluctuation the- 
orems to nonequilibrium systems with chemical reac- 
tions mostly focused on reversible reactions and reac- 
tion networks [U [H, H, [H SH US, HI]. However, 
many reaction-diffusion models are characterized by irre- 
versible reactions, thus yielding a breaking of the usually 
assumed microscopic reversibility. By breaking micro- 
scopic reversibility we mean that if uj(Ci — > Cj) is the 
transition probability from configuration Cj to configu- 
ration Cj, it can happen that w{Ci — > Cj) = even 
though u(Cj — ► C,) > 0. A direct consequence of the 
absence of microscopic reversibility is that many observ- 
ables discussed in the context of fluctuation theorems 
are then ill defined, as in their derivation one explicitly 
uses that for any pat h in configuration space the reversed 
path also exists [H, H|j]. For this reason we focus in 
our study of fluctuations in reaction-diffusion systems on 
an observable that is well defined even in the absence 
of microscopic reversibility. For a system initially in an 
equilibrium steady state this quantity is identical to the 
work observable used in the Jarzynski and Crooks rela- 
tions [Ji, 03. 

It should be noted that diffusion-limited systems with 
effective irreversible reactions can be prepared through a 
fast evacuation of some of the reaction products. This 
makes plausible a possible future verification of the in- 
triguing features that are revealed in our study. 

In the following we discuss, using a numerical ex- 
act method and numerical simulations, steady state and 
transient properties of various reaction-diffusion systems. 
Especially, we study fluctuations in systems that are ini- 
tially in a steady state before being driven away from 
stationarity by varying one of the reaction rates. This 
protocol allows us to measure the probability distribu- 
tion of our observable when going from a steady state 



A, characterized by the value ta of some reaction rate r, 
to another steady state B, characterized by the value tb 
of the same reaction rate. Defining the reversed process 
as changing the reaction rate backwards from tb to r^, 
we can measure also the probability distribution in that 
case and compare the distributions for the forward and 
reversed processes. Even though no exact detailed fluc- 
tuation theorem is observed for the studied quantity in 
the absence of detailed balance, we show that the fluc- 
tuation ratios display intriguing signatures due to the 
specific dynamics of the nonequilibrium systems under 
investigation. 

A brief account of some of our results has been given 
previously 36]. In the present paper we not only give 
a detailed study of the features observed in the ratio of 
the probability distributions for our observable, we also 
extend our investigation to other reaction schemes not 
studied previously. This allows us to gain a better under- 
standing of fluctuations in truly nonequilibrium systems 
driven away from stationarity and to make the discussion 
in [3o| more quantitative. 

Our paper is organized in the following way. In the 
next Section we introduce the different reaction-diffusion 
models and discuss their steady state properties. In Sec- 
tion III we drive these systems out of stationarity through 
time-dependent reaction rates. Using an observable that 
remains well defined even in the absence of microscopic 
reversibility, we study the probability distributions for 
this variable and show that the fluctuation ratios formed 
by the probability distributions computed in the forward 
and the reversed processes display features which can be 
related to the dynamical properties of the nonequilibrium 
system. Finally, Section IV gives our summary as well as 
an outlook on open problems. 



II. THE MODELS AND THEIR STEADY STATE 
PROPERTIES 

We consider one-dimensional lattices made up of L 
sites with periodic boundary conditions where every lat- 
tice site can at most be occupied by one particle. Because 
of the exclusion of multiple occupancy of the lattice sites, 
a total of 2 L configurations exist. Particles are allowed 
to jump to unoccupied nearest neighbor sites with a dif- 
fusion rate D and undergo various creation and annihi- 
lation reactions. We discuss in the following four basic 
reaction schemes, see Table HI and we denote with model 
1, 2, 3, and 4 the four models that result from these reac- 
tion schemes. In all four models we have an annihilation 
process, that takes place with rate A, as well as a cre- 
ation process where a new particle is created with rate h. 
The different models differ by the way creation and anni- 
hilation take place. Let us first discuss the annihilation 
process. In models 1 and 2, two particles on neighboring 
sites undergo a reaction which leads to the destruction of 
one of the particles. This is different in model 3, where 
both particles are destroyed at the same time. Finally, 



in model 4 three neighboring particles are destroyed in 
the annihilation reaction. For the creation process we 
note that whereas in models 2, 3, and 4 new particles are 
spontaneously created at empty sites, in model 1 a new 
particle can only be created at an empty site if one of the 
neighboring sites is already occupied. 



model 1 


model 2 


model 3 


model 4 
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A + A «=» + A 
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A 
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<=> A 
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TABLE I: The different reaction schemes discussed in this 
work. Whereas model 1 is an equilibrium model, in models 
2, 3, and 4 microscopic reversibility is partly or fully broken. 
In the modified models 2', 3' and 4' we allow for reversible 
reactions with rates Ehh and e\X, with < eh < 1 and < 
£a < 1. 

These creation and annihilation processes have been 
chosen in such a way that the models present different de- 
grees of microscopic reversibility. Thus in model 1 all re- 
actions are reversible, and it is easy to see that this model 
is in chemical equilibrium for fixed values of the reaction 
and diffusion rates. Indeed, if both A and h are different 
from zero, the reaction scheme reduces to a unique re- 
versible reaction and detailed balance is fulfilled. Model 
2, on the other hand, does allow for some reactions to be 
irreversible. For example, a new particle can be created 
in the middle of two empty sites, 000 — ► 0A0, with rate 
h, but it is not possible to go back immediately to three 
empty sites as the newly created particle needs a neigh- 
bor for the annihilation process to take place. Finally, 
all reactions are irreversible in models 3 and 4, yielding 
a complete absence of microscopic reversibility. 

We also studied variants of model 2, 3, and 4, called 2', 
3', and 4', where we restore microscopic reversibility by 
allowing the reversed processes to take place with rates 
Ehh and e\X, where < Eh, £\ < L Even though all 
reactions are now reversible, these models do not fulfill 
detailed balance for Eh, £\ < 1 and are therefore still 
nonequilibrium models. 

For all our models the dynamics is described by a 
discrete-time master equation for the probability P(Cj , t) 
that the system is in configuration Ci at time t [37| : 

P(C h t + l)-P(C h t) = Y,W C i -^Ci)P{C h t)- 

j 

u(Ci — »Q)P(Ci,t)] . (1) 

Zia and Schmittmann (38l . |39| pointed out that for this 
type of systems a nonequilibrium steady state is char- 
acterized by both the stationary probability distribution 
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P s {Ci) and the stationary probability currents 

K*(C h C } ) = w(C, — » Ci) P a (C,)-w(C< — > C,)P s (a) 

between two configurations Ci and Cj . 

The stationary probabilities are readily obtained for 
fixed reaction and diffusion rates by setting up the tran- 
sition probability matrix W whose elements are the tran- 
sition rates between different configurations. Of course, 
as we have 2 L configurations for a system of size L, the 
transition probability matrix is a 2 L x 2 L matrix. The 
stationary probabilities are then obtained as the elements 
of the null eigenvector of the Liouville matrix L which re- 
sults when subtracting off the identity matrix from the 
transition probability matrix. For systems that are small 
enough this eigenvalue problem can be solved using stan- 
dard algorithms. For larger system sizes, the stationary 
probabilities can be measured through standard Monte 
Carlo simulations. 

We show in Fig. [1] the stationary probability distribu- 
tions for some of the models and various values of the 
reaction rates. Configurations with the same number of 
particles are grouped together, with the empty configura- 
tion to the left and the fully occupied lattice to the right. 
The first thing to note is that a change of reaction rates 
has a large impact on the stationary probability distri- 
butions. When the creation of new particles takes place 
with a small rate, see the black lines in Fig. [TJ configu- 
rations with only few particles are the most likely. This 
is different when the creation rate is large, as then con- 
figurations with a large number of occupied sites have an 
increasing weight, see the cyan (light gray) lines in Fig. 
[TJ A corresponding behavior is observed when changing 
the rate A. On the other hand, however, a change of 
the value of the diffusion constant D mainly changes the 
distributions quantitatively see Fig. [JJc) and (d). 




configuration C. configuration C. 



FIG. 1: (Color online) The stationary probabilities for (a) 
model 1, (b) model 2, and (c,d) model 3. The common pa- 
rameter is A = 1, whereas D — 1 in (a,b,c) and D = 5 in (d). 
The system size is L = 8. The configurations are grouped by 
the total number of particles in the system, with the empty 
configuration on the left and the fully occupied lattice on the 
right. 



Even though there are visible differences in the sta- 
tionary probability distributions between Fig. [JJJa) and 
(c) , it is not possible to guess from these stationary prob- 
abilities alone whether the system is in equilibrium, as it 
is the case for model I, or whether we are dealing with 
a nonequilibrium system with a fully irreversible reac- 
tion scheme, as for model 3. The stationary probability 
distribution does not allow by itself to characterize un- 
equivocally nonequilibrium steady states. 

Instead of analyzing one-by-one the stationary proba- 
bility currents for various configuration pairs, it is more 
convenient to look at the global quantity [4(| 

K= l**(C<,Q)l ■ (3) 

In Fig. Ufa) we show the dependence of K on the value 
of the creation rate h for fixed values of A and D. Ob- 
viously, the dependence is very different for the different 
models. In the equilibrium model I one does not have 
any non- vanishing stationary probability currents, and K 
is zero for all values of the reaction and diffusion rates, 
as expected. This is different for the nonequilibrium sys- 
tems which are characterized by non-vanishing station- 
ary probability currents. Interestingly, the value of K 
decreases in model 2 for larger h, whereas for models 3 
and 4 it increases as a function of h. In order to under- 
stand this difference in behavior, we recall that for larger 
values of h configurations with a large number of parti- 
cles have an increased stationary probability. As a result, 
free sites will have with high probability occupied neigh- 
boring sites, and the creation process — > A effectively 
equals the process A — > 2A. This is, however, exactly the 
reversed reaction to the annihilation process of model 2, 
which explains why for large h the behavior of model 2 
approaches that of an equilibrium system. For models 3 
and 4, however, all reactions remain irreversible and K 
keeps on growing. 




FIG. 2: (Color online) The total probability current K (a) as 
a function of the creation rate h for models 1, 2, 3, and 4, and 
(b) as a function of e — ex = Eh for models 2', 3', and 4'. In 
all cases A = 1 and D = 1. In (b) the creation rate is h — 2. 
The data are for systems with L = 8 lattice sites. 
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In Fig. ]2jb) we plot K as a function of the parame- 
ter e = E\ = Eh for the modified models 2', 3', and 4' 
for which we allow the reversed reactions to take place. 
As a result, all reactions are now reversible, but for all 
cases we are still out of equilibrium as long as e < 1. In- 
creasing e decreases the distance to equilibrium which is 
finally reached for e = 1 when all reactions and the cor- 
responding reversed reactions are taking place with the 
same rate. 

This discussion of the stationary probability distribu- 
tions and of the stationary probability currents shows 
that it is not possible to characterize in an unambiguous 
way a nonequilibrium system solely through its station- 
ary probability distribution. Much information is con- 
tained in the stationary probability currents which do 
allow to distinguish between the properties of equilib- 
rium, weakly nonequilibrium, and strongly nonequilib- 
rium systems. Therefore these currents allow to quantify 
the distance to equilibrium, making them a useful tool 
for the characterization of nonequilibrium systems. 



III. TRANSIENT FLUCTUATION RELATIONS 

Having discussed the steady-state properties of 
reaction-diffusion models, we now focus on the charac- 
terization of the transient behavior when the systems are 
brought out of stationarity and are then allowed to relax 
to a new steady state. We are realizing this through a 
protocol in which we change one of the reaction rates. 
Experimentally, a change of rates of chemical reactions 
can be achieved by changing the temperature, for exam- 
ple. In our protocol we change one of the rates r from 
an initial value r$ to a final value tm in M equidistant 
steps of length Ar, yielding for the reaction rate the val- 
ues Ti — tq + iAr with i — 0, ■ • • , M. We assume that 
at every step only one reaction or diffusion process takes 
place. 

In the following we discuss mainly numerical exact re- 
sults for small one-dimensional systems. This numerical 
exact approach is rather straightforward and is summa- 
rized in the Appendix. Larger systems can be studied 
along the same lines through numerical simulations, but 
this must be done with some care in order to guarantee 
a sufficient sampling of rare events [4l| . 



A. Observables 

We discuss in the following two observables which dif- 
fer by the fact that for one of the variables microscopic 
reversibility has to be assumed whereas for the other no 
such assumption has to be made. For a system that is mi- 
croscopic reversible, as for example a system that fulfills 
detailed balance, we will discuss the difference between 
these two quantities explicitly. 

In order to define these quantities let us first suppose 
that the system is in a steady state. Starting from a 



configuration Co, the system is in the configuration Cj at 
step i, such that after M steps the system has performed 
the following path in configuration state: X = Co — ► 
Ci — ► • • • — ► Cm-i — > Cm- The probability for this 
path is 



M-l 



P(X) = P S (C ) J] u(Ci 



C 



(4) 



where w(Cj — ► C+i) is the transition probability from 
configuration C to configuration Cj+i. Denoting the re- 
versed path by X = Cm — * Cm-i — > • • ■ — ► Ci — > 
Cp. one then defines for Markovian systems the quantity 



= ln^=ln 

p(x) 



Ps(Cq) 

Ps(C M ) 



M-l ,„ 
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i+1, 



i+1 



Ci 



When the system is driven out of stationarity, we can 
generalize this definition to a time dependent reaction 
rate, yielding 



R = In 



P s (C ,r ) 
Ps{C M ,r M ) 



M-l 

E 

i=0 



111 



LoiCi Cj-|_l 



'i+1 



w(C, 



i+1 



Ci , 7"j ) 



(6) 



where P s (Ci,r i ) is the probability to find the configura- 
tion Ci in the stationary state corresponding to the value 
Ti of the reaction rate r and uj(Ci — > C+i,rj + i) is the 
transition probability from C, to Cj + i at step 

A closer look at the observable R reveals that its def- 
inition requires that if w(Cj — > Cj+i, > than 
w(Cj+i — ► Cj,r,-) also has to be non zero. However, 
in some of our reaction-diffusion models this condition 
is not fulfilled as microscopic reversibility is broken, and 
we can not use R to study them. Hatano and Sasa fl2T ] 
have proposed a different quantity that is closely related 
to R, but that does not assume microscopic reversibility. 
Adapting this quantity for systems driven out of station- 
arity, we can write it in the following way [361 ] 



M-l 
i=0 



(7) 



_P s (Ci,r i+1 ) 
has been called the driving entropy pro- 



The quantity 
duction in [421 ] . 

For a system with microscopic reversibility we can de- 
rive a relation between R and <fi. With the help of the 
probability current 



K*(Ci,C i+ i,ri +1 ) = uj(C i+ i -> Ci,r i+ i)P s (C i+ i,r i+ i 
uj(Ci -> Ci + i,r i+1 )P s (Ci,r i+ i) 

we can write Eq. ((7]) for <j> in the following form: 



(8) 



= R- 
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K*(C 4 ,C 4+1 ,r l+1 ) 



w(C 



i+1 



P s (Cj+i, n+i) lo{C 1+ i 

Ci,r i+1 )~ 
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C u n) 



(9) 
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which reveals that the difference between R and 4> is com- 
posed of terms which have very different physical origins. 
The first term in the In in Eq. ([9]) is due to non- vanishing 
probability currents between different configurations and 
is therefore characteristic for nonequilibrium states. The 
second term is non-trivial only in transient processes as 
it accounts for a shift in the reversed transition probabil- 
ity. This term reduces to the trivial value 1 in case one 
remains in a given steady state, with r,+i = r.; = r$ for 
all i. If this steady state is in addition an equilibrium 
state, the probability currents are all vanishing, and one 
has R = (j> = 0. 

It is easy to show [12, [HI, 52 that for transient pro- 
cesses both quantities fulfill an integral fluctuation the- 
orem: (e~ R ) = 1 and (e~^) = 1, where the average is 
taken over all possible histories when driving the sys- 
tem out of a general steady state. For a system that 
is initially in an equilibrium steady state the relation 
(e~^) = 1 reduces to the Jarzynski relation [J] as then 
(f> = (3(W — AF), where W is the work done on the 
system, AF is the free energy difference between initial 
and final states, and [3 is the inverse temperature. The 
difference Wd — W — AF is the dissipative work. 



B. Probability distributions 

In systems with detailed balance, an exact fluctuation 
relation is obtained when plotting the ratio between the 
probability distributions Pp{(3Wd) and Pr{— @Wd) of the 
dissipative work Wd done on the system in the forward 
and reversed processes [Icj : 



P F {(3W d ) 

p R {-pw d ) 



= e 0W d 



(10) 



Recalling that for systems with detailed balance we have 
the identity <f> — Wd, it is tempting to ask whether for <p 
an exact fluctuation theorem like (|10[) can also be encoun- 
tered for a system initially in a nonequilibrium steady 
state. In fact, this is not the case: the absence of de- 
tailed balance in a nonequilibrium steady state entails 
non-zero probability currents, and no simple relation like 
the relation (fTT)|) exists for cf> in this case. As we shall 
discuss below, the corresponding fluctuation ratios yield 
systematic deviations from the simple behavior encoun- 
tered in systems with detailed balance, these deviations 
containing non-trivial information on the nonequilibrium 
system at hand. 

However, before analyzing these ratios of probability 
distributions, we shall first discuss the probability distri- 
butions themselves. 

Figures EK] show typical examples for the probability 
distributions of R and <f> when changing the creation rate 
from an initial value ho to a final value Km m M steps 
(we only show the case of a varying creation rate h, but 
the following discussion can be made along similar lines 
when changing the value of the annihilation rate A). 
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FIG. 3: (Color online) Probability distributions for the quan- 
tity R when the creation rate is changed in M=6 equidistant 
steps from 0.2 to 1.4 (Pf(R), black curve) or from 1.4 to 0.2 
(Pr(—R), green (gray) curve). The data have been obtained 
for a system with L = 8 sites, with D — 5 and A = 1. (a) 
Model 1, (b) model 2' with e = 0.1, and (c) model 3' with 
£ = 0.1. 



Fig. [3] shows the probability distributions of R for 
three cases that fulfill microscopic reversibility: models 
1,2', and 3'. These different probability distributions are 
not Gaussian but are characterized by a rather irregular 
structure. Their shape depends on the dynamics of the 
different models, expressed by the different reaction 
schemes. It is, however, not straightforward to relate 
specific features of the probability distributions to the 
different reactions. It is important to note that the 
peaks dominating these distributions do not have their 
origin in the noisiness of some numerical data, but are 
real as we are using a numerically exact method. In 
addition, our numerically exact method also allows us 
to circumvent any issues that might appear due to an 
insufficient sampling of rare events. This will be of 
importance in the next section when we discuss the ra- 
tios of the forward and reversed probability distributions. 




FIG. 4: (Color online) The same as in Fig. [3] but now for 
model 3' and two different values of the diffusion rate, (a) 
Pf(R) from the forward process and (b) Pr(—R) from the 
reversed process. 



The probability distributions show a strong depen- 
dence on the system parameters. This is illustrated in 
Fig. 2] where we compare for model 3' the distributions 
obtained for different values of the diffusion rate. When 
we increase the diffusion rate, the general shape of the 
probability distribution changes and, in addition, a large 
number of distinct peaks appear. 
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The probability distributions for <j> differ markedly 
from those for R, see Fig. [5l This was expected as the 
main difference between both quantities are the proba- 
bility currents which are non-zero for a system that is 
out of equilibrium. It is only for the equilibrium model 
1 that the distributions for both quantities are similar. 
Interestingly, the probability distributions for <j> for both 
the forward and reversed processes are characterized by 
the presence of prominent peaks. An increase of the dif- 
fusion constant strongly amplifies these peaks but does 
not change the overall shape of the probability distribu- 
tions. The fact that the heights of the peaks depend on 
the value of the diffusion constant indicates that these 
peaks are related to trajectories in configuration space 
that are dominated by diffusion steps and not by reac- 
tions. In Fig. [5] we verify for model 3 that the main 
contributions to the peaks for a drive of length M = 6 
indeed come from the trajectories where only diffusion 
takes place such that the number of particles is constant 
along these trajectories. The subleading contribution, 
also shown in Fig. O comes from the trajectories where 
a single reaction takes place which changes the number 
of particles in the system. Because the peaks are dom- 
inated by trajectories with pure diffusion, the positions 
of the peaks are the same for the forward and reversed 
processes, the leftmost peak resulting from the diffusion 
of a single particle in the system, whereas the rightmost 
peak is due to the diffusion of a single empty site in the 
system. 

Before closing this section, we remark that in [13] sim- 
ilar peaks have been observed in the probability distri- 
butions of the driving entropy production as well as of 
other related quantities in a model for electron transport 
through a single level quantum dot. 
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FIG. 5: (Color online) Probability distributions for the quan- 
tity cj> when the creation rate h is changed in M=6 steps from 
0.2 to 1.4 (Pf(4>), black curve) or from 1.4 to 0.2 (Pr(-0), 
green (gray) curve). The data have been obtained for a sys- 
tem with L = 8 sites, with D = 5 and A = 1. (a) Model 1, 
(b) model 2, (c) model 3, (d) model 2' with e = 0.1, and (e) 
model 3' with e = 0.1. 




FIG. 6: (Color online) Main contributions to the probability 
distributions for <f> in the forward and reversed processes. The 
black lines show the full probability distributions whereas the 
gray lines show the contributions coming from (a,c) trajecto- 
ries in configuration space with only diffusion steps and no 
reactions and (b,d) from trajectories where exactly one reac- 
tion takes place that changes the number of particles in the 
system. The data are for model 3 with D — 10, ho = 0.2, 
Ah = 1.2, and A = 1. The system size is L — 8 and the 
driving length is M = 6. 



C. Fluctuation ratios 

Having just discussed the probability distributions of 
the quantities R and <f>, we now move on and study the 
fluctuation ratios formed by these probability distribu- 
tions. For a system driven out of an initial equilibrium 
state and fulfilling detailed balance, Crooks has shown 
the exact relation (fTO)) to exist between the probability 
distributions of the dissipative work measured in the for- 
ward and time-reversed processes. This remarkable re- 
sult can be extended to systems that are still reversible 
microscopically but that do not fulfill detailed balance 
any more 16]. As illustrated in Fig. [7] for models 2' 
and 3', the ratios of the probability distributions for R 
show a simple exponential dependence on R. The perfect 
exponential obtained from our data nicely validates our 
numerical exact approach. Obtaining a plot of similar 
quality through Monte Carlo simulations is difficult as 
rare events are then hard to measure. 

Even though in the absence of microscopic reversibility 
R is ill defined, this is different for <p as this quantity ex- 
clusively involves the steady-state probabilities, see Eq. 
([7]) . For an equilibrium system <f> fulfills an exact fluctua- 
tion theorem as it then reduces exactly to the dissipative 
work. As shown in Fig. [5] for model 1, an exponential re- 
lation is indeed obtained for all parameter values as well 
as for different driving processes h{t). 

However, for a system with nonequilibrium steady 
states no exponential detailed fluctuation relation is ex- 
pected for 4> as this quantity does not contain the infor- 
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FIG. 7: (Color online) Fluctuation relation for the observable 
R for model 2' and model 3' for different values of the pa- 
rameter e. The parameters in this calculation are ho = 0.2, 
Ah — 1.2, A = 1, and D = 5. The system size is L — 8 and 
the driving length is M = 6. 




FIG. 8: (Color online) Fluctuation relation for the observable 
<f> for model 1 with (a) different values of D and (b) different 
ways of changing the parameter h(t) with D = 1. The driving 
process usually studied in this paper and which yields the 
data shown in (a) is h(t) ~ t. The parameters used in these 
calculations are ho — 0.2, Ah = 1.2, and A = 1. The system 
size is L — 8 and the driving length is M = 6. 



mation on nonequilibrium currents, see Eq. ©. We show 
in Fig. [5] ratios of the probability distributions of tj> for 
models 2 and 3. For model 2 the deviations from the ex- 
ponential are random and no pronounced dependence on 
system parameters, as for example the diffusion rate D, 
is observed. For model 3, however, a qualitatively differ- 
ent behavior is encountered and systematic deviations in 



the form of oscillations are observed. Similar oscillations 
are also observed for model 4 where three neighboring 
particles are destroyed in the annihilation process. Inter- 
estingly, the amplitudes of these oscillations increase for 
increasing diffusion rates. At first one might think that 
this increase in peak height when increasing D should be 
related to the increase of the peaks in the probability dis- 
tributions themselves, see the discussion in the previous 
section. However, this is too simplistic as an increase of 
peak heights in the probability distributions is also ob- 
served for models 1 and 2 for which we do not observe the 
corresponding behavior in the fluctuation ratios. What 
is different between models 1 and 2 on the one hand and 
models 3 and 4 on the other hand is that for the former 
models any change in the forward and reversed probabil- 
ity distributions is compensated when forming the ratio 
(this compensation is exact for model 1 and approximate 
for model 2), whereas for the latter models this compen- 
sation is only partial, such giving rise to peaks also in the 
fluctuation ratios. 
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FIG. 9: (Color online) Fluctuation ratios for the observable 
(f> for (a) model 2 and (b) model 3 and different values of 
the diffusion constant D. Whereas in model 2 only random 
deviations from a simple exponential behavior are observed, 
systematic deviations show up for model 3. This is highlighted 
in (c) and (d) where we subtract (j> from the logarithm of the 
fluctuation ratio. The light gray lines indicate a simple expo- 
nential dependence. The parameters used in this calculation 
are ho = 0.2, Ah — 1.2, and A = 1. The system size is L — 8 
and the driving length is M — 6. 

Before discussing the origin of this difference, let us 
first have for model 3 a closer look at the peaks in the 
fluctuation ratio. We first note that the positions of these 
peaks are not identical to the positions of the extrema in 
the probability distributions (see for example Fig. 
In Table [Til we compare the positions of the maxima and 
minima in the fluctuation ratio with the peak positions in 
the probability distributions. The observed offset means 
that the peaks in the probability distributions for the for- 
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ward and reversed processes compensate each other when 
forming the ratio, but that the compensation is only par- 
tial away from the peaks. Recalling that the peaks result 
from trajectories in configuration space with only diffu- 
sion steps and that trajectories with reactions make up 
the part between the peaks, we can conclude that re- 
actions are responsible for the peaks in the fluctuation 
ratios. In order to verify this assumption we analyzed 
the contributions to the fluctuation ratio coming from 
the different types of trajectories. We show in Fig. \W\ 
that the observed minima and maxima are indeed mainly 
due to the trajectories with a single reaction process. For 
this we compare the fluctuation ratio with the quantity 
Hf(<I>) /IIr {—(j>) where II(</>) is the probability distribu- 
tion for all trajectories having (a) only diffusion steps or 
(b) exactly one reaction process. Obviously, the peaks in 
the latter ratio coincide with the peaks in the fluctuation 
ratio. 



PD maxima 


FR maxima 


FR minima 


-1.63 


-1.78 


-1.5 


-0.61 


-0.72 


-0.38 


0.60 


0.50 


0.82 


2.02 


1.89 


2.25 


3.64 


3.44 


3.84 



TABLE II: Positions of the maxima in the probability distri- 
butions (PD) and of the maxima and minima in the fluctua- 
tion ratio (FR) for model 3, with D = 5, ho = 0.2, Ah = 1.2, 
and A = 1. The system size is L = 8 and the driving length 
is M = 6. 



ratios also show up in some systems where all reactions 
are reversible. 
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FIG. 11: (Color online) Fluctuation relations for model 3' 
and different values of e. The values of the parameters are 
D = 5, ho = 0.2, Ah = 1.2, and A = 1. The system size here 
is L — 10 and the driving length is M — 10. These data have 
been obtained through Monte Carlo simulations. 




FIG. 10: (Color online) Comparison for model 3 of the fluctu- 
ation ratio (black line) with the ratio IIf(0)/II.r(— (j>) (cyan 
(light gray) line) where Tl(<f>) is the probability distribution 
of <f) for all trajectories with (a) only diffusion steps and (b) 
exactly one reaction process. Note that for trajectories with 
only diffusion few values of <j> can be realized. The common 
parameters are h — 0.2, A = 1., M = 6 and L = 8 and D — 5. 

As a second interesting observation we note that the 
oscillations in the fluctuation ratios are not restricted to 
cases where microscopic reversibility is broken but are 
much more widespread. As is shown in Fig. [11] for model 
3' (the same holds for model 4') peaks in the fluctuation 



In order to understand the origin of these oscillations 
we need to go back to the different reaction schemes 
summarized in Table [T] The configuration space of a 
reaction-diffusion system can be thought to be composed 
of smaller units formed by the configurations with a com- 
mon number N of particles. A diffusion step conserves 
the number of particles, thereby connecting two config- 
urations in the same unit. A passage from one unit to 
another always involves a change of particle number and 
is therefore exclusively due to a reaction process. This 
is sketched in Fig. [T2j Keeping this in mind, a funda- 
mental difference emerges between models 1 and 2 on 
the one hand and models 3 and 4 on the other hand. In 
the former systems every reaction changes the particle 
number by 1, A A" = ±1. In the latter systems, how- 
ever, also larger changes in the particle number happen 
in the annihilation process, with AA" = —2 for model 
3 and AA" = —3 for model 4. As a consequence, loops 
in configuration space that connect a unit with constant 
N with itself and that involve reactions will display an 
asymmetry in the number of creation and annihilation 
processes. Thus for model 3 the smallest loop contains 
two creation processes and one annihilation. This effect 
is still present, even though in a weaker form, when we 
add the backreactions and end up with a microscopically 
reversible model like model 3' with a variable number of 
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particles added or subtracted in the different reactions. 
It is this difference in the number of particles created in 
a creation process or destroyed in an annihilation event 
that yields contributions to the probability distributions 
which are not compensated in the fluctuation ratio. 




FIG. 12: (Color online) Schematic plot of the configuration 
space for models 1, 2, and 3 where configurations are grouped 
by the number of particles in the system N. In a diffusion step 
the system goes from one configuration to another in the same 
unit without changing the particle number. Passages between 
different units are due to reaction processes. For models 1 and 

2 one always has that AN = ±1. This is different for model 

3 (and 4) where different reactions yield different changes in 
the number of particles in the system. 



IV. DISCUSSION 

Characterizing the out-of-equilibrium properties of in- 
teracting many-body systems remains one of the most 
challenging tasks in contemporary physics. The recent 
advent of exact fluctuation and work theorems yielded 
some excitement in the community as it indicated a pos- 
sible way of characterizing large classes of nonequilibrium 
systems. 

In our work we try to characterize diffusion- limited re- 
actions both in their nonequilibrium steady state and in 
the transient state when the systems are driven out of sta- 
tionarity. For systems in their steady state we confirm 
the expectation that probability currents allow to dis- 
tinguish between equilibrium and nonequilibrium steady 
states. In addition, they also allow to define a global 
quantity that quantifies the distance to equilibrium. This 
way of characterizing nonequilibrium steady systems re- 
mains valid even when microscopic reversibility is broken, 
as it is the case for many reaction-diffusion systems. 

The situation is more complicated if one wishes to 
characterize reaction-diffusion systems through fluctua- 
tion and work theorems. If one studies a system for 
which microscopic reversibility is fulfilled, one can de- 
fine a work- like quantity, our quantity R, see eq. ([5]), 
for which exact detailed fluctuation theorems not only 
hold in the steady states but are also valid when the sys- 
tem is driven out of stationarity through time dependent 
reaction rates, fn absence of microscopic reversibility, 
however, R can not be used as it is no longer well de- 



fined. Instead we propose to use the driving ent ropy 
production <f>, see eq. 0, initially introduced in [l2l.l42j]. 
as this quantity exclusively uses stationary probabilities 
and therefore remains well defined even in the absence 
of microscopic reversibility. Whereas the driving entropy 
production always fulfills a global fluctuation theorem 
[H, Ell ! it only fulfills a detailed fluctuation theorem for 
systems with equilibrium steady states. At first look, this 
seems to strongly reduce the usefulness of his quantity 
for the characterization of systems with nonequilibrium 
steady states. However, as we showed in this paper, the 
deviations of the fluctuation ratios for cj> from a simple ex- 
ponential behavior do contain non-trivial information on 
the trajectories in configuration space, fndeed, in cases 
where the change in the number of particles is different 
for different reactions, we observe systematic deviations 
from a simple exponential behavior. These deviations, 
which take the form of peaks superimposed on an expo- 
nential, mainly result from trajectories in configuration 
space where exactly one reaction takes place. 

It is not an easy task to quantitatively relate the peak 
heights and the peak positions to the values of the sys- 
tem parameters. For this, a much more in-depth study is 
needed where all the parameters are varied in a system- 
atic way fill ]. 

Whereas the driving entropy production remains a 
well-defined quantity even in the absence of microscopic 
reversibility, we need to mention that in many cases this 
could be a quantity that is difficult to measure as the 
knowledge of the stationary probability distribution of 
the system is required. As a consequence, the practical 
importance of <p could be restricted, especially for exper- 
imental systems where the stationary probability distri- 
bution is often not easily accessible. 

How general are the results found in this work? Based 
on the reaction schemes discussed in this work and given 
in Table HI we expect the peaks to appear in the fluctu- 
ation ratios for (f> for any reaction-diffusion system that 
allows for a variable number of particles to be created 
or destroyed in the different reactions. This also encom- 
passes more complicated systems with two or more parti- 
cle types. In addition, signatures of the same type should 
also be observed for other system classes with a config- 
uration space topology that is similar to that of the the 
reaction-diffusion systems (i.e., composed by groups of 
configurations that are only connected in a very specific 
way) and with a similar asymmetry in the configuration 
space trajectories. An extension of our work along these 
lines is planed for the future. 
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Appendix 

In order to compute the probability distribution of 
4> 1 see (O, when changing some rate r from its ini- 
tial value ro to the final value tm in M steps, with 
Ti = r a + jj( r M —Ti),i = 0,--- M, we first need to know 
the stationary probability distributions for any value rj. 
This is easily done by determining the null eigenvector 
of the Liouville matrix. We then need to generate all 
possible sequences of configurations (paths in configura- 
tion space) X = Co — > C\ — ► • • • — > Cm-i — ► Cm, 
where only one reaction or diffusion takes place at every 
step. Starting form every possible initial configuration, 
we have to built up a tree structure to all the config- 
urations that can be reached in M steps with non-zero 
probability. This is done recursively by a standard depth- 
first search algorithm that ends when we reach the Mth 
step. We now have to attach a probability to every one 
of these generated paths. For this we are multiplying the 
probability to select the initial configuration Co with the 
product of the M transition probabilities: 

M-l 

P F (X) = P s (Co,ro) J] uj(Ci — > C i+1 ,r i+ i) . (11) 

i=0 

Having now determined every path and its probability, 
we need in addition the values of (f> along these different 
paths, which we obtain through the equation 

M-l 

£(X) - OnP.(Ci,r<+i)-lnP.(a,r<)) , (12) 

i=0 

where P s (C'i,ri) is the stationary probability to find the 
configuration C* at the value n of the rate r. Putting 
everything together, the probability distribution is finally 
obtained through the expression 

= J2 p f ( x ) 5 (?( x ) - A ■ ( 13 ) 

X 

In addition to the just described forward process we also 
study the reversed process where we start in configura- 
tion Cm with the rate tm before changing the reaction 
rate in M steps to its final value tq. The probability 



distribution for this process is then 

p R {4>) = £ p r (x) s (4> (x) - <t) (w) 

x 

with 

M-l 

P R (xj = P s (Cm,tm) Y[ u{C M -i — ► CM-i-i,r M -i-i) ■ 

i=0 

(15) 

In Section III we discuss not only the quantity <f> but 
also the quantity R defined by Eq. ©. For this second 
quantity the procedure is exactly the same, only the cal- 
culation of the values of <j> for the different paths has to 
be replaced by the values of R. 

This numerical exact approach is limited to small sys- 
tem sizes L and few steps M, as the number of paths 
grows exponentially with both L and M, see Fig. IIVI 
For example, for L — 6 the number of paths increases 
from 404 for M = 2 to 8.6 10 s for M = 9. 




timesteps M 



FIG. 13: (Color online) Exponential growth of the calculation 
time in function of the number of steps for model 1 with 
L = 6 sites where the creation rate h was changed between 
ho — 0.2 and hu = 1.4. For this calculation we set A = 1.0 
and considered both vanishing (D = 0) and non-vanishing 
(D = 1) diffusion rates. 
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